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Abstract 

Background: Recent theoretical and empirical work points toward a significant role for sex-chromosome linked 
genes in the evolution of traits that induce reproductive isolation and for traits that evolve under influence of 
sexual selection. Empirical studies including recently diverged (Pleistocene), short-lived avian species pairs with 
short generation times have found that introgression occurs on the autosomes but not on the Z-chromosome. 
Here we study genetic differentiation and gene flow in the long-lived greater spotted eagle (Aquila clanga) and 
lesser spotted eagle (A pomarina), two species with comparatively long generation times. 

Results: Our data suggest that there is a directional bias in migration rates between hybridizing spotted eagles in 
eastern Europe. We find that a model including post divergence gene flow fits our data best for both autosomal 
and Z-chromosome linked loci but, for the Z-chromosome, the rate is reduced in the direction from A. pomarina 
to A. clanga. 

Conclusions: The fact that some introgression still occurs on the Z-chromosome between these species suggests 
that the differentiation process is in a more premature phase in our study system than in previously studied avian 
species pairs and that could be explained by a shorter divergence time and/or a longer average generation time in 
the spotted eagles. The results are in agreement with field observations and provide further insight into the role of 
sex-linked loci for the build-up of barriers to gene flow among diverging populations and species. 



Background 

To get deeper insight into the mechanisms behind popu- 
lation differentiation and speciation, a crucial step is to 
understand the genetic underpinnings of reproductive 
isolation. Earlier efforts have revealed that, given the pro- 
portion of the genome covered by the physical size of 
sex-chromosomes, loci located on these chromosomes 
may contribute appreciably more than expected to 
reduced fitness in hybrids [1]. Most of the evidence 
spring from analyses of species with male heterogamety 
(XY-systems), where advanced genetic or genomic tools 
have been available for some time. It has for example 
been shown that regions on the X-chromosome affect 
hybrid fitness in Drosophila crosses or introgression lines 
[2-4], that several genes that cause hybrid breakdown or 
hybrid sterility in Mus species map to the X-chromosome 
[5-8], and also that the relative size of the X-chromosome 
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compared to the autosomes affects the rate whereby 
reproductive isolation evolves in Drosophila [9]. Recently, 
there has been an accumulation of evidence that points 
to that the Z-chromosome plays a correspondingly 
important role in organisms with female heterogamety 
(ZW-systems), for example birds and butterflies. These 
data are mostly based on other methods than species 
crosses and introgression lines but include observations 
to suggest that Z-linked loci are involved in species 
recognition traits, like coloration [10,11], and female 
mate choice preference [10,12,13], as well as in determin- 
ing hybrid viability and sterility [11,14]. 

The reasons for why sex-chromosome linkage may be 
important for genes involved in reproductive isolation 
are manifold. Sex-linkage enables recessive alleles to be 
expressed in the heterogametic sex, enhancing the 
effects of epistatic interactions involving sex-linked 
alleles. This phenomenon is known as Haldane's rule 
[15,16], and the rule is a commonly applied explanation 
to why hybrids of the heterogametic sex suffer more 
severe fitness reduction in interspecific crosses. Sex- 
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linkage may also facilitate the evolution of sex-specific 
or sexually antagonistic traits and similar to the 'fast X- 
effect' occasionally (but not ubiquitously) observed in 
mammals [17], the Z-chromosome evolves faster than 
the autosomes in birds [18,19]. In female heterogametic 
systems in particular, sex-linkage may enhance the effi- 
cacy of sexual selection since the Z-chromosome is 
inherited directly from father to son [20,21], and sex- 
linked loci evolving under sexual selection might cause 
more rapid advancements in the build-up of reproduc- 
tive isolation than if the loci would have been autosomal 
[22]. Furthermore, although sex-linkage does not imply 
a complete lack of recombination and unless there are 
strong sex-biases in the rate of recombination, it is 
expected that the sex-chromosomes would have a 
reduced recombination rate compared to autosomes of 
similar size since the sex-chromosome only recombines 
in one sex. This could facilitate the diversifying effects 
of reinforcement if both a trait locus and the preference 
for that trait are sex-linked [23], as has been found in 
some species [10,13]. Comparing patterns of genetic dif- 
ferentiation and introgression between genomic regions, 
for example between autosomes and sex-chromosomes, 
is therefore an attractive means to get deeper insight 
into which regions that drive reproductive isolation 
between species of interest. 

The greater spotted eagle Aquila clanga and the lesser 
spotted eagle A. pomarina are two closely related, par- 
tially sympatric, Eurasian birds of prey whose ranges 
overlap in eastern Europe. There is no complete repro- 
ductive barrier between the species and extensive inter- 
breeding, prevalently between A. pomarina males and A. 
clanga females, has been detected across the contact 
zone [24]. The hybridization is introgressive, and 
although gene flow occurs in both directions [24], the 
introgression rate of nuclear markers (AFLPs) has been 
estimated to be about ten times higher than for mtDNA 
which, in agreement with Haldane's rule, possibly 
reflects a lower fitness of hybrid females compared to 
hybrid males [25]. 

Here we use a large subset of previously developed 
gene-based sequences [26-28], to investigate patterns of 
genetic differentiation and to quantify gene flow on the 
autosomes and on the Z-chromosome in a spotted eagle 
hybrid zone. As far as we are aware of this is one of 
very few studies of gene flow in species that are long- 
lived and that have long average generation times. The 
results show biases both between chromosome classes 
and between species and further strengthen the idea 
that sex-linked loci might play an important role in the 
build-up of reproductive isolation between species. 
However, while previous studies on short-lived avian 
species pairs with short generation times have found 
evidence for a very reduced, or even complete absence 



of gene flow on the sex-chromosomes [11,29,30], we 
still find some introgression on the Z-chromosome 
between A. pomarina and A. clanga. 

Methods 

Study species 

The spotted eagles are medium-sized long-lived raptors 
with a generation time of approximately 11 years 
[31,32]. These monogamous birds form sparsely distrib- 
uted solitary pairs and use the same nest over consecu- 
tive years [33]. Both species inhabit mosaic landscapes: 
they nest in forest and hunt in open landscape [33]. 
Although the preferred habitats differ somewhat 
between the species, A. clanga being more closely asso- 
ciated to water, the overlapping distribution ranges 
(Figure 1), partly similar habitat usage and behavioral 
similarities still makes interbreeding possible [34]. 

The numbers of both species have decreased during 
the last century [35], and the decline has been particu- 
larly dramatic in A. clanga, whose vast range across Eur- 
asia is occupied only by few thousand pairs with less 
than thousand pairs breeding in Europe [31]. In con- 
trast, populations of A. pomarina are still dense, and in 
most regions significantly outnumber the sparsely repre- 
sented A. clanga [31,35]. Both species are listed in 
Annex I of the EU Directive on the Conservation of 
Wild Birds (EEC/79/409), as well in the IUCN Red List: 
A. clanga as a globally vulnerable species and A. pomar- 
ina as a species of least conservation concern [36]. 

Sampling and DNA extraction 

We studied 14 unrelated individuals (seven A. clanga, 
seven A. pomarina; three males and four females from 
each species). These were collected in the sympatric 
region in eastern Europe; 13 in Estonia and one in 
Poland (Figure 1). Autosomal markers were analyzed in 
five birds from each of the species, but in the analysis of 
Z-chromosomal markers two more birds were added in 
order to obtain the same number of chromosomes. 
Hence, for both autosomal and Z-chromosomal loci, 10 
chromosomes were analyzed in each of the species 

Blood samples were collected from nestlings or from 
trapped adults. DNA was extracted from blood cells 
using proteinase K treatment followed by a salting [37], 
or a phenol-chloroform purification method [38]. Spe- 
cies-specific morphological characters [39,40] were 
investigated carefully and used for species assignment. 
The final pre-assignment was made, however, using 9 
diagnostic SNPs and 20 microsatellites [41] together 
with a large number of comparative samples, within an 
extensive pan-European hybridization survey [24]. All 
individuals were assigned to either of the species with a 
high probability (> 90%) by two Bayesian model-based 
assignment methods [42,43]. 
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Figure 1 Distribution ranges. Global distribution ranges of the greater spotted eagle and the lesser spotted eagle. The hybrid zone is in the 
area where the species ranges overlap. 



Marker development 

We used a large subset of previously developed gene- 
based sequences [26-28]. First, we tested amplification 
of 122 autosomal and 50 Z-linked primer pairs using a 
single high-quality A. clanga sample, and obtained 87 
and 21 single-band amplification products, respectively, 
in other cases we recorded either several products or no 
DNA-band at all. Strongest single-band PCR products 
were selected for sequencing resulting in high-quality 
sequence data for 36 autosomal and 15 Z-chromosomal 
loci to use in further analysis. All sequences included in 
this study have been submitted to GenBank under 
accession numbers JF521998 - JF522099. 

Amplification was performed in 25 (il containing 25- 
50 ng DNA, 0.25 U AmpliTaq Gold polymerase with 1 
x Amplitaq Gold PCR buffer (Applied Biosystems) or 
Hotstar Taq polymerase with 1 x buffer (Qiagen), 2.5 
mM MgC12, 0.5 uM of each primer and 0.2 mM dNTP. 
The PCR profile included an initial heating at 95°C for 5 
min, followed by 35 cycles of 95°C for 30 s, 60°C to 50° 
C for 30 s and 72°C for 1 min, and a final extension at 
72°C for 10 min. During first subset of cycles (10 or 20), 
an annealing temperature was decreased by 0.5°C or 1°C 
for every cycle, whereas for the remaining cycles 50°C 
was used. 

PCR fragments for sequencing were purified by exo- 
nuclease I and shrimp alkaline phosphatase (USB) treat- 
ment at 37°C for 15 min, followed by denaturation at 
80°C for 15 min. Sequencing was performed by DYE- 
namic ET Terminator or BigDye Terminator sequencing 
reagent premix and MegaBACE 1000 or ABI 3750 auto- 
mated capillary sequencer (Amersham Biosciences) 
according to the manufacturer's recommendations. The 
PCR products were purified in AutoSeq96 plates or 



using an X-terminator purification kit (both by Amer- 
sham Biosciences). 

Data analysis 

All sequences were edited in Sequencher (Gene Codes 
Corp.) and aligned with Clustal W [44] as implemented 
in MEGA4 [45]. Using information of intron-exon 
boundaries from the orthologous genes in chicken, 
sequences were cut to only include the intronic part of 
the gene. All sequences were subsequently purged for 
simple sequence motifs with the help of Sputnik http:// 
espressosoftware.com/sputnik/index.html and individual 
haplotypes were resolved with PHASE v2.1 [46,47]. For 
all loci separately we calculated the population genetic 
summary statistics number of segregating sites (S), 
nucleotide diversity (n), Tajima's D (Tajima 1989), F ST 
and minimum number of recombination events (R) 
using DnaSP [48] (Additional file 1). 

A six-parameter isolation-migration model (IM) 
[49,50] was applied to the data to get estimates of the 
level of gene flow for the different chromosome classes. 
The output of IM contains Maximum Likelihood esti- 
mates and the posterior probability distributions for the 
parameters &i (4N e iU, population mutation rate for 
population 1, N e = effective population size), 0 2 (4N e2 u, 
population mutation rate for population 2), 0 A (4N eA u, 
population mutation rate for the ancestral population), x 
(tu, time since divergence), m 1 {m 1 /[i, migration rate 
from population 2 to population 1 when looking for- 
ward in time) and m 2 (m 2 /u, migration rate from popu- 
lation 1 to population 2 when looking forward in time). 
The model assumes neutrally evolving loci and, comple- 
mentary to the analysis of allele frequency distributions 
(Tajima's D [51], see above), we therefore applied a 
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multi-locus HKA-test [52], as implemented by the soft- 
ware HKA http://genfaculty.rutgers.edu/hey/softwar- 
e#HKA, and a Bayesian method (BAYESFST) [53] to 
investigate if any locus showed evidence for directional 
or balancing selection. Selection tests were applied to 
the autosomal and the Z-linked loci, separately. No 
locus showed evidence for selection at the 1% confi- 
dence level and all loci where therefore assumed to 
evolve neutrally and used in subsequent demographic 
analysis. Inter-specific levels of differentiation (F ST , as 
calculated in DnaSP [48]) and the corresponding 
untransformed p-values for each locus, as calculated in 
BAYESTFST [53], are given in the Additional file 1. A 
second major assumption of IM is that there should be 
no intra-locus recombination [49]. Therefore we applied 
the four-gamete-test (minimum number of recombina- 
tion events > 0) [54] as implemented in DnaSP [48] to 
each locus. Only 6% (3/51) of the loci showed signs of 
recombination after applying the four-gamete-test and 
these were cut so that the longest sequence without evi- 
dence for recombination was analyzed. 

In the IM runs, autosomal and Z-linked loci were ana- 
lyzed separately and each class was run with identical 
settings for three independent runs but with different 
random seeds. Each dataset was run with a wide prior 
parameter range (q (0 = 4N e |i for each population) = 0- 
10, m (m/u for each direction) = 0-50, t (tu) = 0-50) 
in an initial run with a burn-in of 5*10 5 followed by 
5*10 steps and the posterior estimates for the para- 
meters from this initial run was used to narrow the 
priors (q (0 = 4N e u for each population) = 0-2, m (m/u 
for each direction) = 0-25, t (tu) = 0 - 25) in two subse- 
quent, longer analyses. These were run for 5*10 7 and 
1*10 cycles with a burn-in of 1*10 steps, respectively. 
There was good agreement in parameter posterior prob- 
ability distributions and maximum likelihood estimates 
between independent runs and therefore we only report 
the values from the longest run for autosomal and Z- 
chromosome data, respectively (Figure 2, Table 1). The 
mixing of chains was satisfactory and in most runs all 
Effective Sample Size (ESS) values exceeded 2000 and 
for no parameter was the ESS lower than 333 (Table 1). 
All six parameter values were scaled with a per site per 
year mutation rate (u) of 1.4*10" [55] and a generation 
time of 11 years [31,32]. The different inheritance 
modes for autosomes (1) and the Z-chromosome (3/4) 
were accounted for by setting their respective inheri- 
tance scalars in the IM input file. For details about the 
parameter calculations, see the introduction to IM docu- 
mentation http://genfaculty.rutgers.edu/hey/ 
software#IM. 

Finally, we compared the model that included all six 
parameters (® lt 0 2 , &a> t, m x and m 2 ), to a simpler 
demographic model that did not allow for post 



population divergence gene flow (irix = m 2 = 0). These 
analyses were conducted using IMa [56] by running an 
initial M-mode run with the identical settings to the IM 
runs (full 6-parameter model) and sampling 3* 10 s gen- 
ealogies followed by a subsequent L-mode run analyzing 
all possible nested models. The significance of the differ- 
ence between models was evaluated by applying likeli- 
hood ratio tests as implemented in the software. 

Results 

General, levels of polymorphism and signs of selection 

We sequenced 36 autosomal (23.2 kb in total) and 15 Z- 
linked (9.4 kb) genes in 10 individuals from each spe- 
cies. In total, we found 97 single nucleotide polymorph- 
isms (SNPs) in A. clanga and 79 SNPs in A. pomarina 
and 52 of these were shared between the species. Only 
four fixed differences were detected and three of these 
were located on the Z-chromosome. None of the genes 
showed evidence for selection, neither in the HKA test 
(Sum of deviations = 39.5, DF = 100, P = 0.99), in the 
allele frequency distributions (Tajima's D), or in the 
Bayesian analysis after correcting for multiple tests. The 
average nucleotide diversity was higher on the auto- 
somes than on the Z-chromosome in both A. clanga 
(autosomes = (mean) 1.3*10~ 3 ± (SD) 3.8*10~ 4 ; Z-chro- 
mosome = 7.4*10" 4 ± 3.2*10~ 4 ) and A. pomarina (auto- 
somes = 9.7*10 4 ± 3.3*10~ 4 ; Z-chromosome = 7.6*10" 4 ± 
2.8* 10~ 4 ), however, the difference was minor in A. 
pomarina (Wilcoxon's test, W = 305, P = 0.47), and 
only close to significant in A. clanga (W = 358, P = 
0.069). As can be read from the figures above, the aver- 
age autosomal nucleotide diversity was higher in A. 
clanga (1.3*10" ) than in A. pomarina (9.7*10" ), how- 
ever this difference was not significant (W = 767, P = 
0.18), and there was no significant difference between 
the species for the Z-chromosome (7.4*10~ 4 in A. clanga 
and 7.6*10" in A. pomarina, Wilcoxon's test, W = 102, 
P = 0.66). 

Genetic differentiation among species 

The overall F ST between A. clanga and A. pomarina was 
0.30 ± 0.29. The level of differentiation was higher for 
the Z-chromosome (0.37 ± 0.38) than for the autosomes 
(0.27 ± 0.25), but the variance was also higher on the Z- 
chromosome and the difference between chromosome 
classes was not significant (W = 253, P = 0.73). 

Gene flow between species 

We ran IMa analyses for autosomal and Z-chromosome 
data separately and evaluated the significance of a 
model with free and independent gene flow in both 
directions to different nested models with restricted 
gene flow (Table 2). For autosomal data, the model with 
gene flow in both directions was significantly better 
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Figure 2 Posterior probability distributions for the IM model. Posterior probability distributions for the full Isolation-Migration model for 
autosomal and Z-linked loci, respectively. The distribution from the longest (100 Million steps) run is plotted for each chromosomal class. In the 
plots of the effective population size (N e ), red = A clanga, blue = A pomarina and green = ancestral population. In the plots of the migration 
rates, red = migration from A clanga to A pomarina and blue = migration rate from A pomarina to A clanga. 



than all models with restricted gene flow in any or both 
directions (2*log likelihood ratio (LLR) range = 5.01 - 
920.34, df = 1-2, P range = < 1.0*10" 6 - 0.025). However, 
the model with unequal gene flow between species was 
not significantly better than the model with similar gene 
flow in both directions (2*log likelihood ratio = 1.19, df 
= 1, P = 0.27). For Z-chromosome data, the model with 
gene flow in both directions was significantly better 
than a model without any gene flow (2*LLR = 235.39, df 
= 2, P < 1.0*10" 6 ) and a model with no gene flow from 



A. clanga to A. pomarina (2*LLR = 11.2, df = 1, P = 
8.0*10' ) but not better than a model with no gene flow 
from A. pomarina to A. clanga (2*LLR = 0.0048, df = 1, 
P = 0.94). 

The full isolation migration model (IM) revealed 
biases in the patterns of gene flow between species and 
among chromosome classes. For autosomal loci there 
was a marginally higher degree of gene flow from A. 
pomarina to A. clanga (HiPt = 6.8*10" ) than from A. 
clanga to A. pomarina (HiPt = 3.9*10" 6 ) (Figure 2, Table 



Table 1 Summary of posterior values from the IM analysis 





A. da N e 


A pom N e 


Anc. N e 


t 


A da > A. pom 


A. pom > A. da 


» A. pom 


» A. da 


Autosomes 


















HiPt 


8595 


2617 


448 


983 


3.9 


6.8 


0.74 


0.39 


LHPD90 


5862 


1364 


78 


623 


1.6 


2.8 


0.21 


0.083 


HHPD90 


12119 


4721 


52389 


27672 


7.6 


16.3 


2.0 


1.7 


ESS 


8282 


3175 


38954 


333 


5111 


3027 






Z-chromosome 


















HiPt 


3279 


9071 


65 


374 


4.5 


0.034 


0.32 


0.0068 


LHPD90 


1412 


4483 


65 


235 


1.3 


0.0011 


0.040 


0.0011 


HHPD90 


6884 


14688 


43406 


27672 


10.4 


3.4 


1.6 


1.09 


ESS 


9201 


12030 


39997 


2517 


9431 


6606 







Table 1 legend: Summary of posterior values from the IM analysis scaled by a per gene per generation mutation rate of 1.4*10~ 9 . Given is the highest point 
estimate for each parameter (HiPt), the lower (LHPD90) and higher (HHPD90) boundaries for the 90% highest posterior density and the effective sample sizes for 
each parameter (ESS). Only values from the two longest (1 Million steps burn-in + 100 Million steps run time) runs are tabulated. Per generation population 
migration rates into (») each species is also given. A. da = Aquila clanga, A. pom = Aquila pomarina, Anc. = ancestral population, N e = effective population size, > 
= gene flow direction (per gene per generation *10" 5 ), t = time since divergence (1000 years). 
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Table 2 Summary of likelihood ratio test statistics for the 
IMa analysis 

Model log (P) 2LLR df P-value 

Autosomes 



ABCOD (ml =0) 


-2.85 


5.01 


1 


0.025 * 


ABCDO (m2 = 0) 


-171.40 


342.11 


1 


< 1.0*10~ 6 


ABCOO (ml = m2 = 0) 


-460.52 


920.34 


2 


< 1.0* 10~ 6 


ABCDD (ml = m2 > 0) 


-0.94 


1.19 


1 


0.27 


Z-chromosome 










ABCOD (ml =0) 


-6.12 


11.20 


1 


8.0*1 0 4 * 


ABCDO (m2 = 0) 


-0.52 


0.0048 


1 


0.94 


ABCOO (ml = m2 = 0) 


-118.21 


235.39 


2 


< 1 .0*1 0 s 


ABCDD (ml = m2 > 0) 


-0.90 


0.76 


1 


0.37 



Table 2 legend: Summary of likelihood ratio statistics (2LLR), degrees of 
freedom (df) and significance level (P-value) when comparing the full 
likelihood model (ABCDE) to four different nested models (ABCOD (ml = 0), 
ABCDO (m2 = 0), ABCOO (ml = m2 = 0), ABCDD (ml = m2 > 0)), using IMa for 
autosomal and Z-linked loci, respectively. 

1). The proportion of genealogies where this direction 
was inferred was 0.81. For the Z-chromosomal loci, the 
estimated amount of gene flow was similar to the auto- 
somal rate in the direction from A. clanga to A. pomar- 
ina (HiPt = 4.5*10" ), but the rate was severely reduced 
in the direction from A. pomarina to A. clanga (HiPt = 
3.4*10" 8 ) (Figure 2, Table 1). The proportion of genealo- 
gies where the rate from A. clanga to A. pomarina was 
higher than the rate from A. pomarina to A. clanga was 
0.93. These patterns are also evident from the scaled 
population migration rate estimates, which were higher 
for the autosomal loci than for the Z-linked loci into 
both species and drastically reduced into A. clanga for 
the Z-linked loci (Table 1). For most parameters in 
most runs we got good convergence and the posterior 
probability tails were zero within the parameter range. 
However, in some runs, the posterior probability tails 
for ancestral population size and the time since diver- 
gence were still increasing at the highest boundary of 
the parameter interval. 

Discussion 

We analysed 36 autosomal and 15 Z-linked introns in 
population samples of greater spotted eagles and lesser 
spotted eagles from a hybrid zone in eastern Europe. 
The observed levels of genetic diversity were lower for 
Z-linked than for autosomal loci in both species. Since 
the Z-chromosome has a lower population size than 
autosomes, the Z:A population size ratio would be 3:4 
in a population with equal amounts of reproducing 
males and females, which is probably a valid assumption 
in these monogamous species. Hence, the expected 
diversity of the Z-chromosome would be approximately 
75% of the diversity observed for autosomes. This is in 
good agreement with the diversity level observed on the 



A. pomarina Z-chromosome (78% of the autosomal 
diversity level), but the diversity is lower than expected 
for the Z-chromosome in A. clanga (57%). It should be 
noted that the effect of a potential male mutation bias 
would result in an increase of the Z:A diversity ratio. 
Given the limited difference in mutation rate between 
chromosomal classes in birds this should probably only 
have a relatively small effect on the diversity levels 
[55,57], but still suggests that the observed Z-chromo- 
some diversities are at least not higher than expected, 
especially in A. clanga. It has been suggested that selec- 
tion might be more intense on the Z-chromosome than 
on the autosomes and that recurrent selective sweeps 
therefore can cause a reduction in diversity below the 
level expected from differences in population size and 
mutation rates only [58]. Recent data indicate that selec- 
tive sweeps affect the diversity levels only at very limited 
distances from the selected site [e.g. [59,60]], and it 
might be unlikely to expect selective sweeps to occur 
frequently enough to reduce diversity levels over the 
entire chromosome. 

Analogous to the loss of genetic diversity for chromo- 
somes with smaller effective population sizes (see para- 
graph above), genes on the Z-chromosome are expected 
to accumulate allele frequency differences and fixed sub- 
stitutions between diverging populations more rapidly 
than genes on autosomes. We did not find evidence for 
this in the spotted eagles. The level of differentiation 
was on average higher for Z-linked genes than for auto- 
somal genes. However, the variance among genes was 
also higher for Z-linked than for autosomal genes and 
the difference between chromosomal classes was not 
significant. 

Previous work on speciation genetics has pointed out 
that genomic regions harboring genes that affect repro- 
ductive isolation (so called 'genomic islands of specia- 
tion') should experience lower rates of interspecific 
recombination than the genome in general [61,62]. 
These islands may subsequently act as drivers of isola- 
tion through epistatic or physical interactions to other 
loci in the genome resulting in expansion of non-recom- 
bining regions and, at the end, complete reproductive 
isolation [63], although some regions may stay porous 
to post divergence gene flow for significant amounts of 
time [64]. Preliminary attempts using a few loci points 
towards regional variation in degree of interspecific 
recombination among sex-linked genes in Passerina 
buntings [65], and more generally, recurrent analyses 
involving several bird species pairs have revealed that 
the rate of introgression is significantly reduced on the 
Z-chromosome compared to the autosomes. When 
studying multi-locus SNP data in the pied flycatcher 
(Ficedula hypoleuca) and the collared flycatcher (F. albi- 
collis), Saetre et al. [11] found no evidence for 
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introgression on the Z-chromosome despite frequent 
introgression on the autosomes. Similarly, Storchova et 
al. [30] analyzed interspecific migration rates between 
closely related nightingale species (Luscinia luscinia and 
L. megarhynchos) and found that gene flow occurred on 
the autosomes but was completely absent from the Z- 
chromosome. Additionally, Carling et al. [29] found evi- 
dence for autosomal gene flow subsequent to the initial 
divergence of the Lazuli bunting (Passerina amoena) 
and the indigo bunting (P. cyanea) but they could not 
reject a strict allopatric model of divergence when ana- 
lyzing Z-chromosome linked loci. In agreement with 
abovementioned studies our analyses showed that intro- 
gression rates were lower on the Z-chromosome than 
on the autosomes from A. pomarina to A. clanga. In 
fact we could only reject a model without post-diver- 
gence gene flow from A. pomarina to A. clanga for 
autosomal genes, not for Z-linked genes. There was no 
reduction in introgression on the Z-chromosome com- 
pared to the autosomes from A. clanga to A. pomarina 
and we could reject the model without post-divergence 
gene flow, indicating that gene flow still occurs on the 
Z-chromosome in that direction. 

Hence, our data suggest that there has been post- 
divergence gene flow in both directions for the auto- 
somes, but predominantly from A. clanga to A. pomar- 
ina, for Z-chromosome linked genes. There are several 
ways in which these results could be interpreted. In 
agreement with the dominance theory [66,67] of Hal- 
dane's rule [15], it could be that the reduced Z-chromo- 
some introgression from A. pomarina to A. clanga is a 
result of incompatibilities between autosomal alleles spe- 
cific to A. clanga and Z-chromosome alleles specific to 
A. pomarina, and less severely between Z-chromosome 
alleles specific to A. clanga and autosomal alleles speci- 
fic to A. pomarina. A perhaps equally likely and not 
necessarily mutually exclusive explanation is that there 
are sex-biases in hybridization rates between the species. 
It is known from field observations that interbreeding 
usually occurs between A. pomarina males and A. 
clanga females whereas A. clanga males interbreed with 
A. pomarina females less frequently [24]. This is sup- 
ported by mitochondrial data which suggests that 
mtDNA is introgressing into A. pomarina [25], and by 
autosomal data indicating backcrossing mostly to A. 
pomarina [24] . If hybridization only occurs between A, 
clanga females and A. pomarina males and if hybrid 
females have low fitness as data suggest [25], it is possi- 
ble that gene flow on the Z-chromosome would be 
restricted to only occur from A. clanga to A. pomarina. 
This uni-directional gene flow could potentially also 
explain the comparatively high diversity levels observed 
for Z-linked genes in A. pomarina (i.e. introgressed A. 
clanga alleles contributing to higher diversity). 



The divergence time between A. clanga and A. 
pomarina has been estimated to be approximately one 
million years [68]. This is roughly similar to compari- 
sons between other avian species pairs where Z-chro- 
mosome introgression seems to be reduced or 
completely absent [11,29,30,69]. Caution should be 
taken, however, since divergence time estimates that 
apply a molecular clock using few loci might deviate 
significantly from estimates based on likelihood analysis 
applying an isolation migration model [e.g. [70]]. The 
divergence time between A. clanga and A. pomarina is 
based on mitochondrial divergence [68] whereas the 
estimates for both the Luscinia nightingales and the 
Passerina buntings are based on multi-locus data and 
the divergence for the mitochondria is considerably dee- 
per for both of these species [29,30]. A possible expla- 
nation for the low degree of reproductive isolation 
between the spotted eagles compared to previously stu- 
died avian taxa is therefore simply that the divergence 
time is significantly shorter for these species. However, 
the notably longer life span and generation times of 
eagles compared to passerines [31] could also have an 
effect. Hybrid incompatibilities are expected to evolve 
slowly [71,72], and the rate probably decreases with 
increasing life-span and generation times due to mainte- 
nance of ancestral polymorphisms [e.g. [73]]. 

Conclusions 

We studied gene flow, genetic differentiation and 
genetic diversity for 36 autosomal and 15 Z-linked 
genes in a spotted eagle hybrid zone in eastern Europe. 
Our data suggest that introgression occurs in both 
directions but that introgression is reduced on the Z- 
chromosome compared to the autosomes from A. 
pomarina to A. clanga. This is one of few studies ana- 
lyzing long-lived bird species and we show that the 
barriers in the spotted eagles are more permeable to 
gene flow than previously studied avian species pairs. 
The data support an important role for sex-linked loci 
in the build-up of barriers to gene flow and supports a 
model where reproductive barriers evolve in a step- 
wise manner rather than instantly. The observation of 
reduced gene flow on the Z-chromosome together 
with data that indicate a larger variance in level of dif- 
ferentiation between sex-linked than between autoso- 
mal loci makes it tempting to suggest that interspecific 
recombination maintains shared alleles at some sex- 
linked loci while more rapid accumulation of fixed dif- 
ferences occurs at loci with restricted interspecific 
recombination. The latter could constitute potential 
'genomic islands of speciation' and are obvious targets 
for subsequent efforts aiming at identifying loci that 
contribute to reproductive isolation between these 
species. 
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